if !d.name eq 'PS' then begin
    device,xsize=17,ysize=15,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end
default, tt, 14.
@eu_read_bcprofl
@rhoprof
;!p.charsize=2.5
!p.multi=[0,2,1] & !p.charsize=1.0 
!x.margin=[7,0] & !y.margin=[5,3]
the1=reform(the[32,*,0])
;the0=2.0868e6 ; sdr_02
the0=3.4976e6  ; sun_hr_02
dtheta=theta[2]-theta[1]

default,nn,5
ntimes = (size(u))[3]
good = (ntimes-nn)+indgen(nn)-1
print,'doing averages over',nn, ' xaver records ...'

ruu = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
rvv = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
rww = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)

uu=fltarr(l,m)
vv=fltarr(l,m)
ww=fltarr(l,m)
oox=fltarr(l,m)
ooy=fltarr(l,m)
ooz=fltarr(l,m)
pt0=fltarr(l,m)   ;potential temperature
pre=fltarr(l,m)   ;pressure
omega=fltarr(l,m)
ux=fltarr(l,m)
uy=fltarr(l,m)

quu=fltarr(l,m)
qvv=fltarr(l,m)
qww=fltarr(l,m)
qwv=fltarr(l,m)
qwu=fltarr(l,m)
qvu=fltarr(l,m)
rsinth=fltarr(l,m)
r2sinth=fltarr(l,m)
rurms2=fltarr(l,m)
frmc=fltarr(l,m)
frrs=fltarr(l,m)
fthmc=fltarr(l,m)
fthrs=fltarr(l,m)

frt=fltarr(l,m)
ftht=fltarr(l,m)
fx=fltarr(l,m)
fy=fltarr(l,m)

urms2t = ruu+rvv+rww
urms2 = (total(urms2t[good]/float((size(good))[1])))

; reynolds stresses
for k=0, l-1 do begin
  for j=0,m-1 do begin
    rsinth[k,j]=rr*r[k]*sin(theta[j])
    r2sinth[k,j]=(rr*r[k])^2*sin(theta[j])
    uu[k,j] = total(u[j,k,good])/float((size(good))[1])
    vv[k,j] = total(v[j,k,good])/float((size(good))[1])
    ww[k,j] = total(w[j,k,good])/float((size(good))[1])
    oox[k,j] = total(ox[j,k,good])/float((size(good))[1])
    ooy[k,j] = total(oy[j,k,good])/float((size(good))[1])
    ooz[k,j] = total(oz[j,k,good])/float((size(good))[1])
    pt0[k,j] = total(ptemp[j,k,good])/float((size(good))[1])
    pre[k,j] = total(p[j,k,good])/float((size(good))[1])
    quu[k,j]=rho[j,k]*total(u2[j,k,good]- u[j,k,good]^2)/float((size(good))[1])
    qvv[k,j]=rho[j,k]*total(v2[j,k,good]- v[j,k,good]^2)/float((size(good))[1])
    qww[k,j]=rho[j,k]*total(w2[j,k,good]- w[j,k,good]^2)/float((size(good))[1])
    qwv[k,j]=total(rwv[j,k,good]- rho[j,k]*oz[j,k,good]*v[j,k,good])/float((size(good))[1])
    qwu[k,j]=total(rwu[j,k,good]- rho[j,k]*oz[j,k,good]*u[j,k,good])/float((size(good))[1])
    qvu[k,j]=total(rvu[j,k,good]- rho[j,k]*oy[j,k,good]*u[j,k,good])/float((size(good))[1])
    rurms2[k,j]=quu[k,j]+qvv[k,j]+qww[k,j]
    ux[k,j]=ww[k,j]*sin(theta[j])+vv[k,j]*cos(theta[j])
    uy[k,j]=ww[k,j]*cos(theta[j])-vv[k,j]*sin(theta[j])
; r components of the angular momentum flux. all normalized with urms2
    frmc[k,j]=rsinth[k,j]*rho[j,k]*uu[k,j]*ww[k,j]/rurms2[k,j]
    frrs[k,j]=rsinth[k,j]*qwu[k,j]/rurms2[k,j]
; th components of the angular momentum flux. all normalized with urms2
    fthmc[k,j] = rsinth[k,j]*rho[j,k]*(uu[k,j])*vv[k,j]/rurms2[k,j]
    fthrs[k,j] = rsinth[k,j]*qwv[k,j]/rurms2[k,j]
    
    frt[k,j] = (frmc[k,j]+frrs[k,j])
    ftht[k,j] = (fthmc[k,j]+fthrs[k,j])

    fx[k,j]=frt[k,j]*sin(theta[j]) + ftht[k,j]*cos(theta[j])
    fy[k,j]=frt[k,j]*cos(theta[j]) - ftht[k,j]*sin(theta[j])

  endfor
endfor

nlev=20
thg=where(th GE -89 AND th LT 89)

rhov=vv*rh2d*rsinth
rhow=ww*rh2d*r2sinth
stream1=fltarr(l,m)
stream2=fltarr(l,m)
for j=0,m-1 do begin
  for k=0,l-3 do begin
    kk=indgen(k+2)
    stream1[k,j] = int_tabulated(r[kk],rhov[kk,j],/double)
  endfor
endfor
  
for k=0,l-1 do begin
  for j=0,m-3 do begin
    jj=indgen(j+2)
    stream2[k,j] = int_tabulated(theta[jj],rhow[k,jj],/double)
  endfor
endfor

stream=smooth(stream1+stream2,2)

nls=12
lev_stream = 2.*max(abs(stream[*,thg]/rsinth[*,thg]))*indgen(nls)/float(nls-1)$
             -max(abs(stream[*,thg]/rsinth[*,thg]))

;lev_streamn = grange(min((stream[*,thg]/rsinth[*,thg])),0,nls)
;lev_streamp = grange(0,max((stream[*,thg]/rsinth[*,thg])),nls)


;DIF. ROTATION
;omega_0=1.e9/(112.*24.*3600.)                ;112
;lev=grange(omega_0-65,omega_0+70,nlev,/log)  ;112
;title='!6a) !8T!D0!N=112d!6'

;omega_0=1.e9/(56.*24.*3600.)                 ;56
;lev=grange(omega_0-65,omega_0+40,nlev,/log)  ;56
;title='!6b) !8T!D0!N=56d!6'

omega_0=1.e9/(tt*24.*3600.)                 ;28
;lev=grange(omega_0-125,omega_0+65,nlev)       ;28
;lev=grange(300.,470,nlev)       ;28
;title='!6c) !8T!D0!N=28d!6'
title='MHD-02'

;omega_0=1.e9/(14.*24.*3600.)                  ;14
;lev=grange(omega_0-100,omega_0+55,nlev)       ; 14
;title='!6d) !8T!D0!N=14d!6'

print,'Profile of omega computed with T=', tt, 'days'

omega[*,1:m-2]=1e9*(uu[*,1:m-2]/(2.*!pi*rsinth[*,1:m-2]))
omega=omega+omega_0
lev=grange(min(omega),max(omega),nlev)       ;28
print,"Omega_0=",omega_0
print,'levels Omega = ',minmax(omega)
xlength=0.45
xinterval=0.05
xpos1=xinterval
xpos2=xinterval+xlength
cpos=[xpos1,0.04,xpos2,0.94]
bpos=[0.20,0.36,0.22,0.65]
contour,omega[*,thg],x[*,thg],y[*,thg],/fi,nl=64,/iso,pos=cpos,$
	title=title,lev=lev,xtitle='!8x!6',ytitle='!8y!6',xstyle=4,ystyle=4
contour,omega[*,thg],x[*,thg],y[*,thg],nl=21,/iso,/over,lev=lev,color=255
colorbar,range=[min(lev),max(lev)],pos=bpos,/vertical,$
         ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.1,title='!7X!6/2!7p!6 (nHz)'
plots,circle(0,0,0.96),thick=3
plots,circle(0,0,0.61),thick=3
plots,circle(0,0,0.718),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=3,li=0
oplot,[0,0],[0.61,0.96],thick=3,li=0
oplot,[0.4,0.4],[-1.,1.],thick=2,li=1
oplot,[0.6,0.6],[-1.,1.],thick=2,li=1
oplot,[0.8,0.8],[-1.,1.],thick=2,li=1
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
; MERIDIONAL FLOW

lev_vv = 2.*max(abs(vv))*indgen(nlev)/float(nlev-1)-max(abs(vv))
lev = 2.*max(abs(vv))*indgen(nlev)/float(nlev-1)-max(abs(vv))
lev =grange(-5.,5.,nlev)
lev_vv=lev
cpos=[xpos1,0.04,xpos2,0.94]
bpos=[0.70,0.36,0.72,0.65]
jj=thg
contour,smooth(vv[*,jj],2),x[*,jj],y[*,jj],/fi,nl=nlev,/iso,pos=cpos,$
	xtitle='!8x!6',xstyle=4,ystyle=4,lev=lev
contour,smooth(stream[*,thg]/rsinth[*,thg],0),x[*,thg],y[*,thg],/over,lev=lev_stream,$
        max_value=0.,min_value=min(lev_stream),c_thick=2,c_linestyle=2,color=255
contour,smooth(stream[*,thg]/rsinth[*,thg],0),x[*,thg],y[*,thg],/over,lev=lev_stream,$
        max_value=max(lev_stream),min_value=0.,c_linestyle=0,c_thick=3,color=255
colorbar,range=[min(lev_vv),max(lev_vv)],pos=bpos,/vertical,$
         ytickformat='(F6.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.1,title='!8u!D!7h!N!6 (m/s)'
jj=where(th GE -88 AND th LT 88)
;partvelvec,smooth(ux[*,jj],2),smooth(uy[*,jj],2),x[*,jj],y[*,jj],/over,fraction=0.08,length=0.1
plots,circle(0,0,0.96),thick=3
plots,circle(0,0,0.61),thick=3
plots,circle(0,0,0.71),thick=2,li=2
oplot,[0,0],[-0.96,-0.61],thick=3,li=0
oplot,[0,0],[0.61,0.96],thick=3,li=0
print,';;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;'
;
!p.multi=0
save,filename='velocity.sav',omega,vv,ww,x,y,r,th
end
